home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / science / fastmap.zip / TABLE.C < prev    next >
C/C++ Source or Header  |  1993-04-01  |  12KB  |  445 lines
  1. /* Written by Dave Curtis, modified by Gary Williams, HGMP, Harrow, UK */
  2.  
  3. #ifndef __PYRAMID__
  4. #define __PYRAMID__ 0
  5. #endif
  6.  
  7. #include <stdio.h>
  8. #include <math.h>
  9. #if ! __PYRAMID__
  10. #include <string.h>
  11. #include <stdlib.h>  /* prototype for exit() */
  12. #else
  13. char *index();
  14. #define strchr(s,c) index(s,c)
  15. #endif
  16.  
  17. #ifndef max
  18. #define max(x,y) (((x)>(y))?(x):(y))
  19. #endif
  20.  
  21. #define MAXPEDS 100
  22. #define MAXLODS 100
  23.  
  24. #define MIN_RF_TO_FIT 0.009
  25. #define MAX_RF_TO_FIT 0.41
  26. /* for FASTMAP */
  27.  
  28. #ifndef TRUE
  29. #define TRUE 1
  30. #define FALSE 0
  31. #endif
  32.  
  33. int TABLE=TRUE,HOMOG=FALSE,GRAPH=FALSE,BTEST=FALSE,INPUT=FALSE,ATEST=FALSE;
  34.  
  35. #ifndef __ZTC__
  36. #define __ZTC__ 0
  37. #else
  38. unsigned _stack=60000;
  39. #endif
  40.  
  41. #if ! __ZTC__
  42. /* if library does not have strupr() and strstr() */
  43.  
  44. #include <ctype.h>
  45.  
  46. char *strupr(s)
  47. char *s;
  48. {
  49. while (*s)
  50.   {
  51.   if (isascii(*s) && islower(*s))
  52.       *s=toupper(*s);
  53.   ++s;
  54.   }
  55. }
  56.  
  57. char *strstr(s,t)
  58. char *s, *t;
  59. {
  60. while ((s=strchr(s,*t))!=0)
  61.     if (!strncmp(s,t,strlen(t))) return s;
  62.     else ++s;
  63. return 0;
  64. }
  65.  
  66. #endif /* __ZTC__ */
  67.  
  68. error(e,s1,s2)
  69. int e;
  70. char *s1,*s2;
  71. {
  72. fprintf(stderr,"%s%s\n\n",s1,s2);
  73. exit(e);
  74. return 0;
  75. }
  76.  
  77. usage()
  78. {
  79. #ifdef UNIX
  80. printf("\nRun on FILENAME.RES (OUTFILE.DAT) file from MLINK or LINKMAP\n\n\
  81. Usage:\n\
  82. TABLE [-g -h -b -i -a] filename.res\n\n\
  83. -g  -> make .grp file for EASIPLOT\n\
  84. -h  -> make .hom file for HOMOG\n\
  85. -b  -> make .bte file for BTEST\n\
  86. -i  -> make .inp file for FASTMAP\n\n\
  87. -a  -> add A-test results to table and graph files\n\n\
  88. (files for HOMOG and BTEST will only be created if results are\n\
  89. available for multiple pedigrees)\n\n");
  90. #else
  91. printf("\nRun on FILENAME.RES (OUTFILE.DAT) file from MLINK or LINKMAP\n\n\
  92. Usage:\n\
  93. TABLE filename.res [/g /h /b /i /a]\n\n\
  94. /g  -> make .GRP file for EASIPLOT\n\
  95. /h  -> make .HOM file for HOMOG\n\
  96. /b  -> make .BTE file for BTEST\n\
  97. /i  -> make .INP file for FASTMAP\n\n\
  98. /a  -> add A-test results to table and graph files\n\n\
  99. (files for HOMOG and BTEST will only be created if results are\n\
  100. available for multiple pedigrees)\n\n");
  101. #endif
  102. return 0;
  103. }
  104.  
  105. int LINKMAP=0,MLINK=1,basepoints=0;
  106.  
  107. float kosambi(theta)
  108. double theta;
  109. {
  110. return (theta>0.49 ? 70.0 : 25*log((1+2*theta)/(1-2*theta)));
  111. }
  112.  
  113. #define MAXTHETAS 10
  114.  
  115. main(argc,argv)
  116. int argc;
  117. char *argv[];
  118. {
  119. FILE *fi,*fo,*fg;
  120. float logli0[MAXPEDS],theta[MAXLODS][MAXTHETAS],*lod[MAXPEDS],
  121. totlod[MAXLODS],fdum,f,
  122. theta0[MAXTHETAS],almax[MAXLODS],amax[MAXLODS];
  123. float distance[MAXLODS];
  124. char filename[81],s[81],*ptr,title[81],famid[MAXPEDS+1][31];
  125. int i,nval,nfams,idum,j,lorder[MAXTHETAS],testlocus,locusnum,rrvary,nloci;
  126. if (argc<2) {usage();error(1,"","");}
  127. strcpy(filename,argv[1]);
  128. while (--argc>1)
  129.   {
  130.   strupr(argv[argc]);
  131.   if (!strcmp(argv[argc],"-G") || !strcmp(argv[argc],"/G")) GRAPH=TRUE;
  132.   else if (!strcmp(argv[argc],"-H") || !strcmp(argv[argc],"/H")) HOMOG=TRUE;
  133.   else if (!strcmp(argv[argc],"-B") || !strcmp(argv[argc],"/B")) BTEST=TRUE;
  134.   else if (!strcmp(argv[argc],"-I") || !strcmp(argv[argc],"/I")) INPUT=TRUE;
  135.   else if (!strcmp(argv[argc],"-A") || !strcmp(argv[argc],"/A")) ATEST=TRUE;
  136.   else {usage();error(2,"Bad switch: ",argv[argc]);}
  137.   }
  138. #ifndef UNIX
  139. if (!strchr(filename,'.'))
  140.    strcat(filename,".res");
  141. #endif
  142. fi=fopen(filename,"r");
  143. if (!fi) error(2,"Could not open ",filename);
  144. #ifdef UNIX
  145. if (!strchr(filename,'.'))
  146.    strcat(filename,".");
  147. #endif
  148. nfams=nval=0;
  149. *title='\0';
  150. strcpy(title,"\n");
  151. while (fgets(s,81,fi))
  152.    {
  153.    strupr(s);
  154.    if ((ptr=strstr(s,"MLINK"))!=0 && strchr(s,':'))
  155.       strcpy(title,ptr);
  156.    if ((ptr=strstr(s,"LINKMAP"))!=0)
  157.       {
  158.       if (strchr(s,':'))
  159.          strcpy(title,ptr);
  160.       MLINK=0;
  161.       LINKMAP=1;
  162.       }
  163.    if (LINKMAP && strstr(s,"TEST LOCUS"))
  164.       sscanf(strchr(s,':'),": %d",&testlocus);
  165.       /* linkmap */
  166.    if (MLINK && strstr(s,"VARY"))
  167.       {
  168.       sscanf(strchr(s,':'),": %d",&rrvary);
  169.       --rrvary;
  170.       }
  171.    if (strstr(s,"ORDER"))
  172.      {
  173.      nloci=sscanf(strchr(s,':'),": %d %d %d %d %d %d %d %d %d %d",
  174.        &lorder[0],&lorder[1],&lorder[2],&lorder[3],&lorder[4],&lorder[5],
  175.        &lorder[6],&lorder[7],&lorder[8],&lorder[9]);
  176.      for (locusnum=0;locusnum<MAXTHETAS;++locusnum)
  177.          if (lorder[locusnum]==testlocus) break;
  178.      break;
  179.      }
  180.    }
  181.  
  182. while (fgets(s,81,fi))
  183.    {
  184.    strupr(s);
  185.    if (LINKMAP && (ptr=strstr(s,"LOCUS ORDER"))!=0)
  186.      {
  187.      sscanf(strchr(s,':'),": %d %d %d %d %d %d %d %d %d %d",
  188.        &lorder[0],&lorder[1],&lorder[2],&lorder[3],&lorder[4],&lorder[5],
  189.        &lorder[6],&lorder[7],&lorder[8],&lorder[9]);
  190.      for (locusnum=0;locusnum<MAXTHETAS;++locusnum)
  191.        if (lorder[locusnum]==testlocus) break;
  192.      }
  193.    if ((ptr=strstr(s,"THETAS"))!=0)
  194.    {
  195.    if (basepoints==0)
  196.      {
  197.      ++basepoints;
  198.      sscanf(ptr,"THETAS %f %f %f %f %f %f %f %f %f %f ",
  199.        &theta0[0],&theta0[1],&theta0[2],&theta0[3],&theta0[4],&theta0[5],
  200.        &theta0[6],&theta0[7],&theta0[8],&theta0[9]);
  201.      while(fgets(s,81,fi),strupr(s),(ptr=strstr(s,"LIKE"))==0) ;
  202.      fgets(s,81,fi);
  203.      while (fgets(s,81,fi),
  204.            sscanf(s," %29s %f %f",famid[nfams],&fdum,&logli0[nfams])==3)
  205.               if (++nfams>MAXPEDS) error(4,"Too many pedigrees - increase MAXPEDS","");
  206.      if (!nfams)
  207.        {
  208.        while (fgets(s,81,fi),strupr(s),
  209.          !strstr(s,"TOTALS")) ;
  210.        sscanf(s,"TOTALS %f %f ",&fdum,&logli0[0]);
  211.        }
  212.      }
  213.    else
  214.      {
  215.      if (basepoints==1)
  216.        {
  217.        ++basepoints;
  218.        for (i=0;i<(nfams?nfams:1);++i)
  219.          {
  220.          lod[i]=calloc(MAXLODS,sizeof(float));
  221.          if (!lod[i]) error(5,"Out of memory - try decreasing MAXLODS","");
  222.          }
  223.        }
  224.      sscanf(ptr,"THETAS %f %f %f %f %f %f %f %f %f %f ",
  225.        &theta[nval][0],&theta[nval][1],&theta[nval][2],&theta[nval][3],
  226.        &theta[nval][4],&theta[nval][5],&theta[nval][6],&theta[nval][7],
  227.        &theta[nval][8],&theta[nval][9]);
  228.      if (MLINK) distance[nval]=kosambi(theta[nval][0]);
  229.      else if (locusnum==0) distance[nval]= -kosambi(theta[nval][0]);
  230.      else
  231.        for (distance[nval]=0,i=0;i<locusnum;++i)
  232.           distance[nval]+=kosambi(theta[nval][i]);
  233.      totlod[nval]=0.0;
  234.      while(fgets(s,81,fi),strupr(s),(ptr=strstr(s,"LIKE"))==0) ;
  235.      fgets(s,81,fi);
  236.      for (i=0;fgets(s,81,fi),
  237.            sscanf(s," %29s %f %f",famid[i],&fdum,&f)==3;++i)
  238.               {
  239.               lod[i][nval]=f-logli0[i];
  240.               totlod[nval]+=lod[i][nval];
  241.               }
  242.      if (!nfams)
  243.        {
  244.        while (fgets(s,81,fi),strupr(s),
  245.          !strstr(s,"TOTALS"))  ;
  246.        sscanf(s,"TOTALS %f %f ",&fdum,&lod[0][nval]);
  247.        lod[0][nval]-=logli0[0];
  248.        }
  249.      ++nval;
  250.      }
  251.    }
  252.    }
  253. fclose(fi);
  254.  
  255. if (MLINK)
  256.   {
  257.   for (i=0;i<nfams;++i)
  258.     lod[i][nval]=0;
  259.   totlod[nval]=0;
  260.   theta[nval][rrvary]=0.5;
  261.   ++nval;
  262.   }
  263.  
  264. for (i=0;i<nval;++i)
  265.   {
  266.   if (totlod[i]<-99.0) totlod[i]= -99.0;
  267.   for (j=0;j<nfams;++j)
  268.     if (lod[j][i]<-99.0) lod[j][i]= -99.0;
  269.   }
  270.  
  271. if (ATEST && nfams>1)
  272.   {
  273.   for (i=0;i<nval;++i)
  274.     {
  275.     float alpha,alod;
  276.     amax[i]=0.0;
  277.     almax[i]=0.0;
  278.     if (i<nval-1)
  279.       for (alpha=0.05;alpha<1.01;alpha+=0.05)
  280.       {
  281.       alod=0.0;
  282.       for (j=0;j<nfams;++j)
  283.          alod+=log10(alpha*pow(10.0,lod[j][i])+1-alpha);
  284.       if (alod>almax[i])
  285.         {
  286.         almax[i]=alod;
  287.         amax[i]=alpha;
  288.         }
  289.       }
  290.     }
  291.   }
  292.  
  293. if (TABLE)
  294. {
  295. strcpy(strchr(filename,'.'),".tab");
  296. fo=fopen(filename,"w");
  297. if (!fo) error(2,"Could not open ",filename);
  298. if (*title) fprintf(fo,"%s\n",title);
  299. if (LINKMAP)
  300.  {
  301.  for (j=0;j<nloci-1;++j)
  302.   {
  303.   fprintf(fo,"theta  ");
  304.   for (i=0;i<nval-1;++i)
  305.    fprintf(fo,"%8.3f%s",theta[i][j],(i==nval-2)?"\n":"");
  306.   }
  307.  fprintf(fo,"cM     ");
  308.  for (i=0;i<nval-1;++i)
  309.    fprintf(fo,"%8.3f%s",distance[i],(i==nval-2)?"\n":"");
  310.  }
  311. else 
  312.   {
  313.   fprintf(fo,"theta  ");
  314.   for (i=0;i<nval-1;++i)
  315.    fprintf(fo,"%8.3f%s",theta[i][rrvary],(i==nval-2)?"\n":"");
  316.   fprintf(fo,"cM     ");
  317.   for (i=0;i<nval-1;++i)
  318.    fprintf(fo,"%8.3f%s",kosambi(theta[i][rrvary]),(i==nval-2)?"\n":"");
  319.   }
  320. for (j=0;j<nfams;++j)
  321.   {
  322.   fprintf(fo,"%6.6s ",famid[j]);
  323.   for (i=0;i<nval-1;++i)
  324.     fprintf(fo,"%8.3f%s",lod[j][i],
  325.       (i==nval-2)?(nfams==1?"\n\n\n":"\n"):"");
  326.   }
  327. if (nfams>1) for (i=0;i<nval-1;++i)
  328.     fprintf(fo,"%s%8.3f%s",i?"":"total  ",
  329.       totlod[i],(i==nval-2)?"\n\n\n":"");
  330. if (ATEST && nfams>1)
  331.   {
  332.   for (i=0;i<nval-1;++i)
  333.     fprintf(fo,"%s%8.3f%s",i?"":"alpha  ",
  334.       amax[i],(i==nval-2)?"\n":"");
  335.   for (i=0;i<nval-1;++i)
  336.     fprintf(fo,"%s%8.3f%s",i?"":"A-lod  ",
  337.       almax[i],(i==nval-2)?"\n\n\n":"");
  338.   }
  339. fclose(fo);
  340. }
  341.  
  342. if (INPUT && !LINKMAP)
  343. {
  344. strcpy(strchr(filename,'.'),".inp");
  345. fo=fopen(filename,"w");
  346. if (!fo) error(2,"Could not open ",filename);
  347. for (j=0;j<nfams;++j)
  348.   {
  349. /*
  350.      fprintf(fo,"%6.3f 0.01 %6.3f 0.2 %6.3f 0.4\n",
  351.          lod[j][2],lod[j][5],lod[j][7]);
  352. */
  353.     {
  354.     for (i=0;i<nval;++i)
  355.        if (theta[i][rrvary]>=MIN_RF_TO_FIT && theta[i][rrvary]<=MAX_RF_TO_FIT)
  356.          fprintf(fo,"%6.3f %7.4f ",theta[i][rrvary],max(lod[j][i],-99.0));
  357.     fprintf(fo,"\n");
  358.     }
  359.   }
  360. fclose(fo);
  361. }
  362.  
  363. if (GRAPH)
  364. {
  365. strcpy(strchr(filename,'.'),".grp");
  366. fg=fopen(filename,"w");
  367. if (!fg) error(2,"Could not open ",filename);
  368. fprintf(fg,"PARMS:L\n");
  369. fprintf(fg,"TITLE:%s",title);
  370. fprintf(fg,"TITLEV:lod\nTITLEG:%s\nTITLEG:%s\nTITLEG:%s",
  371.    LINKMAP?"cM":"theta",
  372.    LINKMAP?"cM":"theta\nTITLEG:cM",
  373.    nfams==1?"lod score":"total");
  374. if (nfams>1) for (j=0;j<nfams;++j)
  375.    fprintf(fg,"\nTITLEG:ped %s",famid[j]);
  376. if (ATEST && nfams>1) fprintf(fg,"\nTITLEG:alpha\nTITLEG:atotal");
  377. if (LINKMAP)
  378.   for (i= -1;i<nval;++i)
  379.    {
  380.    fprintf(fg,"\n%f,%f",i==-1?-70:distance[i],
  381.      (i==-1)?0.0:totlod[i]);
  382.    if (nfams>1) for (j=0;j<nfams;++j)
  383.       fprintf(fg,",%f",(i==-1)?0.0:lod[j][i]);
  384.    if (ATEST && nfams>1)
  385.       fprintf(fg,",%f,%f",(i==-1)?0.0:amax[i],(i==-1)?0.0:almax[i]);
  386.    }
  387. else
  388.   for (i=0;i<nval;++i)
  389.    {
  390.    fprintf(fg,"\n%f,%f,%f",theta[i][rrvary],
  391.      kosambi(theta[i][rrvary]),
  392.      totlod[i]);
  393.    if (nfams>1) for (j=0;j<nfams;++j)
  394.       fprintf(fg,",%f",lod[j][i]);
  395.    if (ATEST && nfams>1)
  396.       fprintf(fg,",%f,%f",amax[i],almax[i]);
  397.    }
  398. fprintf(fg,"\nDATATYPE:XY1,%d,-1\n",LINKMAP?2:3);
  399. if (nfams>1) for (j=1;j<=nfams;++j)
  400.   fprintf(fg,"DATATYPE:XY1,%d,-%d\n",j+2+(LINKMAP?0:1),j+1+(LINKMAP?0:1));
  401. fclose(fg);
  402. }
  403.  
  404. if (HOMOG && nfams>1)
  405. {
  406. strcpy(strchr(filename,'.'),".hom");
  407. fo=fopen(filename,"w");
  408. if (!fo) error(2,"Could not open ",filename);
  409. fprintf(fo,"%4d\n1\n",nval-1);
  410. for (i=0;i<nval-1;++i)
  411.   fprintf(fo,"%6.3f%s",MLINK?theta[i][rrvary]:distance[i],(i==nval-2)?"\n":" ");
  412. for (j=0;j<nfams;++j)
  413.   {
  414.   fprintf(fo,"1\n");
  415.   for (i=0;i<nval-1;++i)
  416.     fprintf(fo,"%6.3f%s",lod[j][i],(i==nval-2)?"\n":" ");
  417.   }
  418. fprintf(fo,"0\n0\n");
  419. fclose(fo);
  420. }
  421.  
  422. if (MLINK && BTEST && nfams>1)
  423. {
  424. strcpy(strchr(filename,'.'),".bte");
  425. fo=fopen(filename,"w");
  426. if (!fo) error(2,"Could not open ",filename);
  427. fprintf(fo,"50\n%d\n",nfams);
  428. fprintf(fo,"%d\n",theta[0][rrvary]?nval-1:nval-2);
  429. for (i=theta[0][rrvary]?0:1;i<nval;++i)
  430.   fprintf(fo,"%6.3f%s",theta[i][rrvary],(i==nval-1)?"\n":"");
  431. for (j=0;j<=nfams;++j)
  432.   for (i=theta[0][rrvary]?0:1;i<nval;++i)
  433.     if (lod[j][i]<-9.9) lod[j][i]= -9.9;
  434. for (j=0;j<nfams;++j)
  435.   {
  436.   for (i=theta[0][rrvary]?0:1;i<nval;++i)
  437.     fprintf(fo,"%6.3f%s",lod[j][i],(i==nval-1)?"\n":"");
  438.   }
  439. fclose(fo);
  440. }
  441. return 0;
  442. }
  443.  
  444.  
  445.